Unconventional magnetism via optical pumping of interacting spin systems 
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We consider strongly interacting systems of effective spins, subject to dissipative spin-flip processes 
associated with optical pumping. We predict the existence of novel magnetic phases in the steady- 
state of this system, which emerge due to the competition between coherent and dissipative processes. 
Specifically, for strongly anisotropic spin-spin interactions, we find ferromagnetic, antiferromagnetic, 
spin-density-wave, and staggered- XY steady states, which are separated by nonequilibrium phase 
transitions meeting at a Lifshitz point. These transitions are accompanied by quantum correlations, 
resulting in spin squeezing. Experimental implementations in ultracold atoms and trapped ions are 
discussed. 



Exotic magnetic states play a central role in the physics 
of quantum many-body systems, and have been explored 
in a wide variety of strongly correlated materials (ij. 
Realizing and exploring magnetic states has recently 
emerged as a central goal in ultracold atomic physics 
[2, yj. Due to highly controllable and tunable interac- 
tions, ensembles of ultracold neutral atoms and ions may 
provide a unique laboratory to study exotic quantum 
magnetism [2|-l8| Among the main obstacles are relatively 
small energy scales associated with magnetic ordering 
(e.g., the superexchange scale in the Hubbard model), 
requiring cooling atomic systems down to very low tem- 
peratures [2| and the slow timescales involved in spin 
thermalization [9l4ll|. Furthermore, ultracold atoms are 
fundamentally open, driven quantum systems far away 
from their absolute thermal equilibrium. This motivates 
the exploration of spin dynamics in the presence of driv- 
ing and dissipation 12|-|29f . 

Recently a number of schemes involving dissipation to 
create magnetic phases have been proposed. These typ- 
ically use engineered reservoirs involving coupling multi- 
ple lattice sites 12h14|. At the same time, one expects 
single-site dissipation such as spontaneous decay to be 
detrimental to realizing interesting magnetic states, re- 
sulting e.g. in unwanted decoherence. In this Letter, 
we demonstrate that optical pumping and spontaneous 
decay can instead enrich the phase diagram, resulting in 
new phases and phase transitions that do not exist in con- 
ventional equilibrium systems. Significantly, these novel 
states can be observed under conditions when realization 
of conventional, equilibrium states is difficult. 

The key idea of this work can be understood by con- 
sidering the anisotropic spin- 1/2 Heisenberg model (i.e., 
the XYZ model), which is governed by the Hamiltonian 



^ = ^ E (-^-^ + Jv« + J^«)^ (1) 



d-dimensional cubic lattice with nearest-neighbor inter- 
actions. In the presence of conventional optical pumping, 
this Hamiltonian is augmented with a dissipative process 
that flips the spins down at some rate 7 (i.e., it corre- 
sponds to the jump operator a~ on every site, where 

The steady state of this open many-body system is 
easy to understand in the case of isotropic spin-spin 
interactions, namely the XXZ model (with either fer- 
romagnetic or antiferromagnetic couplings). For this, 
the Hamiltonian can be rewritten in the form H = 
(l/2d) E[2J.(^+ff,7 + ^™^+) + J.«]. This Hamil- 
tonian conserves the total number of spins in the | f) 
state, and therefore does nothing to counteract the spon- 
taneous decay. Thus, the steady state is a trivial dark 
state with all spins polarized, | || • • • |) (|| • • • | | , so the 
XXZ model never experiences a phase transition in the 
presence of dissipation, regardless of Jx and Jz- 

However, new types of magnetic order emerge 
for strongly anisotropic couplings. The crucial 

role of anisotropy can be understood as follows. 
Each spin experiences an effective magnetic field 
{Jx{<J^), Jy{ay), Jz{a^)), which depends on the direction 
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where a!!ll,a^,a^ are the Pauli matrices for an effective 
spin n. We assume that the spins are localized on a 
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FIG. 1. Mean-field phase diagrams for the dissipative XYZ 
model with (a) Jz/"/ — 1 and (b) Jz ~ 0, showing the different 
phases: paramagnetic (PM), ferromagnetic (FM), antiferro- 
magnetic (AFM), spin-density- wave (SDW), and staggered- 
XY (sXY). The white arrow points to a Lifshitz point. 



of its neighbors [Fig. [U^a)]. It precesses about this ef- 
fective field and also decays towards | |). In order for 
the spin to point away from | I) in steady state, its 
precession must be strong enough to counteract the de- 
cay. In the isotropic case, the spin is always parallel 
to the magnetic field, so there is no precession at all. 
On the other hand, when the couplings are sufhciently 
anisotropic (e.g., Jx ~ ~Jy)^ the spin is roughly per- 
pendicular to the magnetic field, so the precession is 
strong enough to point the spin away from | \.) [Fig.[2lja)]. 
This is in sharp contrast with thermal equilibrium state, 
in which the spin tries to align with the magnetic field 
rather than precess about it. 

This competition between precessional and dissipa- 
tive dynamics gives rise to a remarkable phase diagram 
(Fig. [T]), including ferromagnetic and antiferromagnetic 
phases as well as spin-density-wave and staggered- XY 
phases that do not exist in equilibrium. The spin-density- 
wave, paramagnetic, and ferromagnetic phases meet at 
multicritical Lifshitz points, at which the period of the 
spin-density wave diverges [30| ; such Lifshitz points have 
been seen in equilibrium magnets with long-range inter- 
actions [3l|, |32], but generally do not exist in nearest- 
neighbor spin models. In addition, we find that a contin- 
uous symmetry emerges for certain couplings; the spon- 
taneous breaking of this symmetry leads to a phase we 
call the staggered- XY phase. Finally, we find that quan- 
tum correlations (as measured by spin squeezing) persist 
near the phase transitions. 

The model described here can be implemented in sys- 
tems of trapped ions or systems of ultracold atoms with 
anisotropic superexchange or dipolar interactions. The 
spin states | f) f^^d | \) correspond to two electronic 
states of the ion or atom. In the case of ions, the spin- 
spin interaction is obtained through virtual transitions 
involving motional sidebands [J, |33J, |3J| . In the case of 
ultracold atoms, the spin-spin interaction is obtained us- 
ing a two-photon resonance that excites and de-excites 
atoms in pairs [35|, as explained in the Supplementary 
Material, or using superexchange interactions in p-band 
optical lattices [36|. In all cases, dissipation can be con- 
troUably introduced using optical pumping. 

Model. We now turn to detailed analysis of the phe- 
nomena outlined above. The dynamics of the many-body 
system are given by a master equation for the density 
matrix p: 



P 



-i[H,p\+lY, 



^n P^n " 2 (f^rV« P + P^u^n ) 



(2) 



Equation ([2]) has a unique steady state solution [37[ , and 
we are interested in whether the steady state exhibits 
a phase transition as the parameters J^;, Jj,, Jz change. 
Note that the decay is independent for each spin, in con- 
trast with the Dicke model 28, bsl. Furthermore, the 



spins are not in equilibrium with the environmental bath. 
Thus, in contrast with the spin-boson model 39|,|40|, the 
steady state is not the joint ground state of the system 
and environment. 

The master equation has a Zi symmetry (cr^ , ct^ — > 
— (7^,— cr|(), which will be spontaneously broken by the 
magnetic phases, as explained below. In practice, there 
may also be dephasing noise, leading to dissipative terms 
in Eq. ([2]) like a^pa^; since the Z2 symmetry is unaf- 
fected by these terms, the phase transitions we describe 
are robust to dephasing, although the phase boundaries 
are shifted. 

Mean-field theory. We begin by solving for the steady 
states of the model Eq. ^ at the level of mean-field 
theory. We allow the mean field to vary on each site 
to account for spatially inhomogeneous states [20]. The 
mean-field equations, which are simply nonlinear Bloch 
equations, are: 

m 

^ = -7(«> + l) 

+ \Y.^UO{0 - JyiOiOl (3) 
m 

where the sum over m is taken over nearest neighbors 
of n. (A related model with only dephasing noise was 
studied in Ref. 4l|, |42| . Another related model with an 
external field and nonlinear damping was studied using 
the Landau-Lifshitz-Gilbert equation [43|, |4J|.) 

Clearly, there is always a fixed-point solution, (cr^) — 
(a^) ~ 0, (crfj) = —1, in which all the spins are pointing 
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FIG. 2. (a) Bloch-sphere plot, showing mean-field values of 
(ct) (solid red arrow) and effective magnetic field (dashed blue 
arrow) for Jxh = '^Jy/l = fj J^ = 0- The vectors are 
normalized to unit length, (b), (c) show the sXY phase in the 
xy plane of the Bloch sphere, (b) shows one possible stable 
configuration. Black and red arrows correspond to sublattices 
A and B. (c) shows that the A sublattice (black solid arrow) 
generates a magnetic field (black dashed arrow) that the B 
sublattice (red solid arrow) precesses around. Similarly, the 
B sublattice generates a magnetic field (red dashed arrow) 
that the A sublattice precesses around. The angle 6 can take 
any value. 



down. We call this the paramagnetic (PM) phase, since it 
does not break the Z2 symmetry of Eq. ([2]) • We now con- 
sider the linear stability of the PM phase as a function of 
JxjJy^Jz '45!|. We consider d-dimensional perturbations 
with wave vector k = (fci, ^2, • ■ • , ka) where kg, = I-k j ag 
and at is an integer. We find that the PM phase is un- 
stable to perturbations of wave vector k when 



:^|:cosfc,-j.^<-g.(4) 



This condition is satisfied only when the couplings are 
sufficiently anisotropic. 

When the PM phase is unstable, the system ends up 
in a time-independent steady state with {a^, {alf^) ^ 0, 
so it breaks the Z2 symmetry of the master equation. 
There are four types of ordered phases: (i) Spatially uni- 
form state, which we call the ferromagnetic (FM) phase, 
resulting from instability of the PM phase to A:^ = for 
all (.. (ii) Spatially modulated state with a period of 
two lattice sites in all directions, i.e., the system divides 
into two sublattices. We call this the antiferromagnetic 
(AFM) phase, and it results from instability to A:^ = tt 
for all L (iii) Spatially modulated state with a period 
greater than two lattice sites in at least one direction, 
which we call the spin-density-wave (SDW) phase. This 
results from instability to all other kg. (iv) When J^ = 0, 
there is also a staggered-XY (sXY) phase, resulting from 
instability to all fcf , which is discussed below. The phase 
diagram is shown in Fig. [Ija)-(b). The transitions from 
the PM phase are continuous, while the FM-AFM tran- 
sition is discontinuous. 

We note two unusual features of this phase diagram. 
First, along the boundary between the PM and SDW 
phases, the k at which the instability of the PM occurs 
approaches 0, meaning that the period of the SDW di- 
verges [Fig. EJa)]. This line culminates in a multicriti- 
cal Lifshitz point between the PM, FM, and SDW 
phases. Lifshitz points occur in magnetic models with 
competing interactions 3lll32|. but are not found in equi- 
librium nearest-neighbor magnets: thus, their existence 
in nearest-neighbor magnets out of equilibrium indicates 
that nonequilibrium phase diagrams can be qualitatively 
richer than those in equilibrium. Lifshitz points show 
enhanced fluctuation effects relative to conventional crit- 
ical points |30| , and hence offer a rich venue for studying 
quantum fluctuations away from equilibrium. 

The second distinctive feature of the phase diagram 
is that the ordered phase breaks a continuous symme- 
try when Jz = 0. In this case, the system divides into 
two sublattices like in the AFM phase. However, the 
angle between the two sublattices can take any value. 
In the specific case of J^; = —Jy, the spins on the A 
and B sublattices are at angles 9 and —9 relative to the 
X = y line on the Bloch sphere [Fig. [^fb)]. Any value of 
9 corresponds to a stable configuration, since the sublat- 




FIG. 3. (a) Unstable wave vector k along the lower boundary 
of the PM phase in Fig. [ija). A Lifshitz point occurs at 
Jx/'y = 1.32. For convenience, only one-dimensional wave 
vectors are shown, (b) Squeezing parameter (" , calculated 
in the Gaussian approximation for J^ = 0. The sXY phase 
has been whited-out, since the Gaussian approximation is not 
valid there. 



tices remain perpendicular to each other's magnetic field 
[Fig. mjc)]. Upon ordering, this continuous C/(l) sym- 
metry between the sublattice spin orientations is sponta- 
neously broken, leading to a phase we call the staggered- 
XY (sXY) phase. This phase has vortex-like topological 
defects around which the relative orientation between A- 
and B-sublattice spins rotates by 27r. 

Comparison with equilibrium. It is instructive to con- 
trast the above results with the equilibrium case (for 
d > 1). The equilibrium ground state of Eq. ([1]) is or- 
dered for any J^, Jy, Jz |46[. The magnetization axis is 
determined by the strongest of the coupling constants, 
and the sign of that coupling determines whether the or- 
dering is ferromagnetic or antiferromagnetic. Evidently, 
the nonequilibrium phase diagram exhibits qualitatively 
different behavior from this equilibrium case. The qual- 
itative differences between equilibrium and nonequilib- 
rium remain even in the limit 7 — > 0, although the steady 
state takes an increasingly long time to reach. 

Fluctuation effects. We now turn from mean-field 
theory to an analysis of fiuctuations. Such an analy- 
sis was recently performed for driven polariton conden- 
sates [47| and suggests that the static critical proper- 
ties (i.e., renormalization-group fixed points) of a driven 
Markovian system are related to finite-temperature equi- 
librium critical properties. This would indicate that the 
dissipative XYZ model discussed here undergoes true 
phase transitions in two or more dimensions. 

We estimate fiuctuation effects and squeezing in the 
Gaussian approximation by mapping the spins to hard- 
core bosons dil: o-+ -> b'l,a^ -> 26jj6„ - 1. This 
gives a reliable approximation in the PM phase, where 
(ffi) ~ — 1- To Gaussian order (which includes relaxing 
the hardcore constraint), the resulting Hamiltonian is 



H 



1 
2d 



iJ^ + Jy)T.(''^^nbn + bmbi) 

{inn) 

(J. - Jy) Y. (f'l^n + bn.br.) ~ 4dJ. ^ 6j,6„ 



(5) 
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FIG. 4. Correlation function (a^a^) for ID chain of 16 spins, 
from simulating the master equation, (a) Jz/7 — 1, showing 
remnant of FM for Jx/^ — 2, Jy = (blue circles, solid line); 
remnant of AFM for Jx/j = —2, Jy = (green triangles, 
dashed line); remnant of SDW for Jx/^ = 4, Jy/-y — 2 (red 
squares, dash-dotted line). The period of the SDW matches 
the mean-field prediction (5.3 sites), (b) Jx/j ~ —Jy/"/ = 
l, Jz/"/ ~ 0, showing remnant of sXY phase. 



and the dissipative terms in the master equation are 
"fJ2ni^ripbl^ - ^{b'lbnP + pt'lbn)]- We now use standard 
Keldysh path-integral techniques [48| to compute the re- 
laxation rate, (cr^), and the squeezing. We summarize 
the results here and provide details in the Supplemental 
Material. 

(1) Relaxation rate. The rate at which the steady 
state is approached can be read off from the poles of 
the retarded Green's function. For notational simplic- 
ity, we assume d = 1 here. In the Gaussian approxi- 
mation, the lowest pole has complex frequency —ij/2 ± 
2y/{JxC0sk — Jz){JyCOsk — Jz). A continuous phase 
transition occurs when the frequency of this pole ap- 
proaches zero; this precisely recovers Eq. (|4]). 

(2) Below-threshold fluctuations. Near the transition, 
one expects to find nonanalytic behavior in the number 
of up spins, X^n ('''«)■ -'^O'" "^2 — 0, this scales as {(7z) "^ 
(7^ + IQJxJy)'^'^^'^^^'^. The divergence for d = \ renders 
the Gaussian approximation inconsistent, and is related, 
as we shall show in a future work, to the absence of a 
phase transition in one dimension (consistent with the 
poIariton-BEC case (47|). 

(3) Squeezing. We find that spin squeezing, a measure 
of quantum correlations, persists near the transition. It 
can be calculated using the definition of squeezing for 
bosons 0: C' = 1 + 2((6^6) - |(6)|2) - 2|(62) - (6)2|. For 
the case of Jz ~ 0, as the phase boundary is approached, 
C^ — )■ 2 in ths thermodynamic limit for the fc = 0, tt 
modes, signaling the presence of quantum correlations 
[Fig-lib)]. 

Comparison with numerics. We have also simulated 
the E g. ^ in ID using the method of quantum trajec- 
tories [50j . Although there is presumably no phase tran- 
sition in ID, the numerical results already show quali- 
tative features predicted by mean-field theory. For ex- 
ample, when mean-field theory predicts FM, the cor- 
relation {(J^'^n) is positive for all distances [Fig. lUJa)]. 
When there should be AFM, the correlation alternates 



sign. When there should be SDW, the correlation varies 
with a wavelength that matches the mean-field value. 
When there should be sXY, {a'^a'^, (o'm^n)' (^m'^'K) 
are all for odd distances and positive for even distances 
[Fig. mjb)]. Furthermore, the gap of the Liouvillian ap- 
proaches at the boundary of the FM phase, consistent 
with the Gaussian approximation (see Supplemental Ma- 
terial). 

Experimental realization. The dissipative XYZ model 
can be implemented experimentally using trapped ions. 
One can use ^^^Yb"*" and let | J,) and | f) correspond 
to ^S'l/alF = 0,TOi. = 0) and ^1^3/21^^ = 2,mF = 0). 
In the presence of laser beams judiciously detuned from 
certain motional sidebands, the ions interact via Eq. ([T]) 
[J, |33|, |34[ . Jx.Jy, Jz can be on the order of 1-5 kHz, and 
their magnitudes and signs can be varied by changing 
the laser detunings [J|. By admixing a small component 
(10""') of ^P3/2 using an off-resonant laser, one broadens 
the linewidth of | f) to 2 kHz. (To make this a closed 
cycle, additional lasers optically pump back into | J,) on a 
much faster timescale.) Thus, the parameter space shown 
in Fig.[l]is experimentally achievable. This setup can im- 
plement an arbitrary lattice topology for a large number 
of ions [5l| . 

A variety of other realizations of the XYZ model are 
also possible. One approach is to use ultracold atoms 
coupled via dipole-dipole interactions. The XYZ Hamil- 
tonian is implemented by driving a two-photon reso- 
nance, so that atoms are excited and de-excited in pairs, 
as explained in the Supplemental Material. This scheme 
can be realized using Rydberg-dressed atoms [S^, Ryd- 
berg atoms [351, ISSJ, Isj], or dipolar atoms or molecules 
[55|. We show explicitly in the Supplemental Mate- 
rial that, for Rydberg-dressed atoms, the parameters 
needed for the phase transitions (Fig. [T]) are experimen- 
tally achievable. Finally, one can adapt a recent proposal 
for realizing XYZ_ models via superexchange in p-band 
optical lattices [363 to include dissipation, by optically 
pumping the atoms into the Px orbital via an intermedi- 
ate excited orbital (e.g., dx^^y^) that does not decay into 
the s band. 

Conclusion. In summary, we have computed the phase 
diagram of anisotropic spin models subject to sponta- 
neous decay, and shown that these models exhibit phases 
(SDW and sXY) and phase transitions (Lifshitz point) 
that are not found in similar equilibrium models. The 
qualitative differences can be traced to the fact that in 
equilibrium, spins align with the magnetic field, whereas 
away from equilibrium, they precess about it. We find 
that quantum correlations, as measured by squeezing, 
persist near the dissipative transitions. This work paves 
the way for future explorations of critical behavior and 
nonequilibrium fiuctuations near the phase transitions we 
have identified. A particularly intriguing question is how 
frustrated interactions (due to a triangular lattice) affect 
the AFM and sXY phases. 
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Supplemental material: 
Unconventional magnetism via optical pumping of interacting spin systems 

Tony E. Lee, Sarang Gopalakrishnan, and Mikhail D. Lukin 

In this supplement, we first discuss a level scheme for experimentally realizing the XYZ Hamiltonian using dipole- 
dipole interactions. Then we elaborate on the calculation of excitation gaps and fluctuations presented in the main 
text. Finally, we show numerical results for the Liouvillian gap. 



Level scheme for anisotropic spin Hamiltonians 



We outline a general scheme for realizing the anisotropic Hamiltonian J^afj^crf-^ + JyC^^cry + J^cr^'^n starting with 
atoms with Ising interactions, which can originate from dipole-dipole interactions. (Our scheme is similar to that of 
Ref. 35|, but we use two-photon transitions instead of four-photon transitions.) This desired Hamiltonian can be 



rewritten in terms of raising and lowering operators as ^//(c^cr„ -|- cr„a+) + Jsq{(T^a^ + cr^a„ ) -|- JzCr^a^, where 

Jff = J. 

to anisotropic, 



Jy and Jsq 



Jx — Jy The first term corresponds to isotropic, flip-flop interactions; the second corresponds 
"squeezing" interactions that raise or lower a pair of neighboring spins at once. Our primary interest is 
in realizing the regime where Jx ~ —Jy, so that the flip-flop interaction is much weaker than the squeezing interaction. 





FIG. 5. Level structure and laser coupling scheme for realizing anisotropic XYZ spin models. We assume two-level atoms, with 
a transition frequency cj, and an interaction of strength V in the z basis. The levels are coupled using four lasers f^i (red), ^2 
(dark blue), fl^ (orange) and fl^ (light blue), such that the pairs (r2i,r23) and {fl2,Q4,) are on two-photon resonance. We take 
the detunings A to be much smaller than V. The left panel shows processes that involve raising or lowering two spins at a 
time; the right panel shows flip-flop processes. The processes denoted with gray arrows and marked with an X are suppressed 
because the associated energy denominators are large. 

Consider two-level atoms, each with ground state \g) and excited state |e). We explain our scheme assuming that 
the dipole-dipole interaction exists only between the excited states (like with Rydberg atoms), but the scheme also 
works when there are interactions between ground states or between ground and excited states. The bare atomic 
Hamiltonian for two atoms takes the form 



H = uj{a^+a^) + Va'i''a'^ 



(6) 



where uj is the transition frequency from g to e, af^ = |e)(e|m, and V is the level shift due to the interaction between 
the excited states. We now couple the atoms using four lasers ili, J72, ila, J74 as shown in Fig. [5l such that pairs of 
lasers (r2i,J73) and (512,5^4) are on two-photon resonance, and their detunings from the intermediate atomic levels 
are A. Using the assumption that A ^ V, we ignore intermediate states with energy denominators ~ V. Then the 
allowed transitions are those shown in Fig. [Sj from second-order perturbation theory, 

'2A "^ 2A ~ 2A 



Jff 



2A 



(7) 



Note that in the hmit oi A -^ V, Jff and Jgq are independent of V. It follows, in particular, that if we choose all 
four fli to have the same magnitude, and arrange the relative phases so that r^ifia = —fl2^4^ then the flip-flop term 
cancel out completely, but the squeezing terms from the two pairs of lasers add constructively: Jff — 0, Jsq = 2fll/A 
and thus J^ = ^Jy = Vl\/ A. To obtain the J^crf^crfj term in the Hamiltonian, one detunes the lasers from two-photon 
resonance. 

The Ising Hamiltonian Eg. (|^ and thus the two-photon resonance scheme can be realized in a number of different 
settings: Rydberg atoms [35|, |53J, |5J] , Rydberg-dressed atoms 52|, or dipolar atoms or molecules [55|- In the Rydberg 
case, the interactions are naturally Ising in character. Using dipolar molecules, one can realize Ising interactions by 
using the TO = ±1 states in the J — \ rotational manifold [? ]. These levels are separated from the m — state 
by an electric field- induced splitting; within the to. = ±1 manifold, flip-flop interactions are prevented by angular 
momentum conservation, so that the dipolar interaction is purely Ising. 

In each of the cases mentioned above, one can tunably engineer dissipation by dressing one of the spin states off- 
resonantly with a short-lived excited state. For concreteness, we provide example numbers for the Rydberg-dressed 
case. Let | i) and | f) correspond to the F = \ and F = 2 hyperfine levels of the 551/2 ground state of ^'^Rb. By 
dressing | f) with a small component (10^^) of the 50S Rydberg state, the dipole-dipole interaction at a distance 
of 500 nm is 44 MHz [56|. By doing the two-photon scheme with $7^ = 1 MHz and A = 10 MHz, one obtains 
Jx = —Jy — 100 kHz. Then by dressing | f) with a small component of the 5P3/2 excited state, the effective linewidth 
of I t) can be controllably set to, for example, 100 kHz. 



Keldysh calculations in the Gaussian approximation 



In what follows we briefly describe how excitation gaps, fluctuations, and squeezing can be computed using the 
Keldysh technique 28|, |48| . We shall analyze the following master equation: 



i[H,p]+jJ2 



b,pbl~l{pblb, + blb,p) 



(9) 



For notational convenience, we choose dimension d = 1. Then 



H = ^^ T.(blb, + b^b]) + ^V^ E(^l^' + b^b,) - 2 J. Y: bh 



(10) 



iv) 



(u> 



In what follows we define f = J^ + Jy and A = Jx — Jy Using standard identities, one can rewrite the master 
equation as the following Keldysh action, in Fourier space: 



S 



duodkl ^ a\[uj + 2J^ - /cos(fc)][6* (w, fc)6a(w, fc) +b*^{-uj,-k)ba{-uj.,-k)] 



a=± 



-A cos{k)[b*jLu, k)bl{-Lu, -k) + ba{uj, k)ba{-Lo, -k)] 



-17 



-«7 



b+{u;, k)b*_{u}, k) - -[b+{Lj, k)bX{uj, k) + 6_(a;, k)b*_{uj, k)] 



b+{-cj, -k)b*_i-u;, -k) - -[fo+(-a;, -fe)6;(-a;, ~k) + b_{-u}, -fc)61(-w, -k)] 



(11) 

(12) 
(13) 



where ba{uj, k) is now a complex field. If we define the fields bd = (b+ + b-)/\/2 and bq = (6+ — b-)/\/2 (where all 
other arguments are assumed to be the same), then this action can be written in the compact 4x4 form: 



D^'^ 
S = I dujdk{b*i{uj,k) bci{-uj,~k) b*{uj,k) 6,(-w, -fc) ) ( ^^ ^^ 



in terms of the inverse retarded Green's function 

ik + 
-Acosfc 



jf _ I Lu — f COS k + 2Jz + ij/2 



— Acosfc 
-uj — f cos k + 2Jz — J7/2 



/ bci{oj,k) \ 

blii-^.-k) 

bq{uj,k) 

\b;i-u:,-k) J 



(14) 



(15) 



and the inverse Keldysh Green's function 



^K _ I il 



^=(o;)- (i«) 



Relaxation rate 



In order to find the relaxation rate (i.e., the eigenvalue gap of the Liouvillian) we use the fact that the characteristic 
frequencies in the system are given by the zeros of the determinant of D^ . These frequencies correspond to the poles 
of the retarded Green's function [48J , and are easily computed to be 



w = -i-f/2 ± V(7cosAr^2XF^^^cos2fc. (17) 

When the w vanishes, there is a second-order phase transition. This occurs when A^ cos^ k — (/cosfc — ^J^)"^ = 7^/4, 
which, when re-expressed in terms of Jx and Jy, gives the result shown in the main text. 

Expectation values and squeezing 

In order to compute expectation values and correlations in the steady state, it is necessary to compute the so-called 
Keldysh Green's function Gk 28[ |48| , which is given by the matrix product 



Gk = -{D")-^D^{D^^)-\ (18) 

We consider the case oi Jz — 0. Transforming to the time domain and setting t = 0, we find that 



/2(&t(fc)6(fc))+l 2(6t(fc)6t(_fc)) 
iL.K[K)-\ 2{b{k)b{-k)) 2(6t(fc)6(fc))+l 



1 / 7-^+4/^cos^fc -2Acosfc(2/cosA: + i7) 



^2 ^ 4( j2 _ ^2) cos2 k \ -2A cos fc(2/ cos k ~ i7) Y + 4r cos^ ^ 



(19) 



From this, it is straightforward to arrive at the expression for (6^6) given in the main text. Moreover, using the 
definition of squeezing in Ref. [49J, one can compute the squeezing as plotted in the main text. 

Liouvillian gap 

Below are plots of the Liouvillian gap, calculated numerically. The Liouvillian gap is defined as the real part of 
the eigenvalue of the master equation with largest nonzero real part. It corresponds to the rate of approach to the 
steady state. The gap was calculated using exact diagonalization for a ID chain of 6 spins with nearest-neighbor 
interactions. The gap become small at the boundary of the paramagnetic phase (compare with the phase diagram 
in the main text), indicating a critical slowing down. In particular, for the case of J2/7 = 1, gap clearly shows the 
transition from the paramagnetic phase to the ferromagnetic and antiferromagnetic phases. The gap also becomes 
small at the transition from the paramagnetic phase to the spin-density-wave phase (the oblong shape is probably 
due to the small size of the system, since only certain wave vectors are allowed). 

For an infinite system in two or more dimensions, we expect the gap to vanish along the boundary of the paramag- 
netic phase; the fact that the gap does not vanish here is due to the finite size and the low dimensionality. Note that 
the gap is 1/2 along the J^ = Jy line, because the Hamiltonian contains only flip-flop terms {(y1^a~ -I- o'~(7+) that do 
not counteract the overall decay. 
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FIG. 6. Gap of the Liouvillian for (a) Jzj^ = 1 and (b) Jz ~ 0. 



